/*
 * EquilibriumCondition.h
 *
 *  Created on: 8 Aug 2011
 *      Author: Allan
 */

#ifndef EQUILIBRIUMCONDITION_H_
#define EQUILIBRIUMCONDITION_H_

// C++ includes
#include <string>
using namespace std;

// Eigen includes
#include <Eigen/Core>
using namespace Eigen;

// GeoReact includes
#include "../Utils.h"

// GeoReact forward declarations
class EquilibriumSolver;
class Multiphase;
class ReactionSystem;

typedef std::function<const double(double T, double P, const VectorXd& n, const VectorXd& a)> 
EquilibriumConditionFunction;

class EquilibriumCondition
{
public:
	EquilibriumCondition();
	
	EquilibriumCondition(const Multiphase& multiphase, const ReactionSystem& equilibriumReactions);
	
	void
	SetTotalComposition(const string& species, double moles);
	
	void
	SetComposition(const string& species, double moles);
	
	void
	SetActivity(const string& species, double activity);
	
	void
	SetAcidity(double pH);
	
	void
	SetPartialPressure(const string& gas, double Pg);
	
	void
	SetChargeBalance();
	
	friend class EquilibriumSolver;
	
private:
	/// The pointer to the multiphase object which holds all the species and phases
	const Multiphase* multiphase;
	
	/// The pointer to the reaction system object which holds all equilibrium-controlled reactions
	const ReactionSystem* equilibriumReactions;
	
	/// The number of primary (Nj) and equilibrium (Ne) species
	uint Nj, Ne;
	
	/// The indexes of the primary (jIndexes) and equilibrium (eIndexes) species in the multiphase object
	vector<Index> jIndexes, eIndexes;
	
	/// The alpha and beta row that represents mathematically this equilibrium condition
	VectorXd alphaRow, betaRow;
	
	/// The function h that represents mathematically the right-hand side of this equilibrium condition
	EquilibriumConditionFunction h;
	
private:
	const uint
	LocalIndexPrimarySpecies(const string& species) const;
	
	const uint
	LocalIndexEquilibriumSpecies(const string& species) const;
	
	void
	SetCompositionAuxiliar(const string& species, const EquilibriumConditionFunction& h);
};

inline const EquilibriumConditionFunction 
Constant(double val)
{
	return [=](double T, double P, const VectorXd& n, const VectorXd& a) { return val; };
}

#endif /* EQUILIBRIUMCONDITION_H_ */
